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Abstract 

Additive noise is known to produce counter-intuitive behaviors in nonlinear dynamical systems. 
Previously, it was shown that systems with a deterministic limit cycle can display bistable switching 
between metastable states in the presence of asymmetric additive white noise. Here, we systematically 
analyze the dynamics of this bistable behavior and show how the vector field away from the limit 
cycle influences the rate and directionality of the bistable switching. Using stochastic phase reduction 
methods, we identify mechanisms underlying different rates of switching and predict when the system 
will rotate in the opposite direction of the deterministic limit cycle. Thus, this work presents an 
alternative mechanism for generating a range of bistable switch-like behaviors that have been observed 
in a number of physical systems. 
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1 Introduction 

Countless physical systems display dynamics that are modeled by limit cycle oscillators [13,16,23]. As 
many of these systems are subjected to stochastic fluctuations, there has been a growing interest in 
understanding how noise affects the dynamics of limit cycle oscillators [4,10,11,21,24], When the noise 
is weak, formal reduction methods can be performed to reduce the dimensionality of the system to a 
so-called phase equation [21,24], In this case, it is known that the magnitude of the added noise shifts 
the mean frequency of oscillations away from the natural frequency of the limit cycle. On the other 
hand, when the noise is large, the trajectories of the system bear no resemblance to the deterministic 
limit cycle behavior. 

The case of moderate noise applied to a limit cycle oscillator has received less attention. In this case, 
the vector field around the limit cycle has a strong influence on the dynamics of the stochastic system, 
as the trajectory is constantly being perturbed off the limit cycle. In previous work [15], we showed that 
an oscillator subject to moderate noise can display numerous interesting and counter-intuitive behaviors. 
In some cases, the addition of noise can act to completely eliminate oscillations, and even cause the 
trajectories to rotate in the opposite direction of the deterministic limit cycle. Perhaps more interesting 
is the emergence of bistable switching behavior which occurs in noisy non-rotationally symmetric systems. 

Bistable switching behavior has been observed in variety of physical systems, including gene networks 
[9,14], neural systems [5,7], and neural growth processes [2] . In the simplest example, bistable switching 
arises when noise is added to a system that has deterministic bistability, i.e., Kramers’ problem [8] . In this 
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case, the stochastic process will randomly switch between the two deterministic fixed points (left panel 
in Figure 1). The resulting stationary density of the process will have two distinct peaks centered around 
the deterministic stable fixed points (right panel in Figure 1). In [15], we demonstrated that similar 
switching dynamics arise when noise is added to systems with a deterministic limit cycle. Importantly, 
in the examples studied in [15], the limit cycle always exists, and thus the points where the stochastic 
process spent most of its time represent metastable points that do not exist in the deterministic system. 
Through the analysis of two simple planar limit cycle systems, we were able to identify precisely when 
this bistable switching will occur. However, the underlying dynamics of the bistable behavior in these 
systems can be quite variable (see Figure 3). Here, we systematically analyze the bistable behavior in 
these simple systems, and identify the mechanisms leading to different types of switching dynamics. 


Stationary Density 



t 


Figure 1: Bistable switching. The left panel shows a stochastic process w(t) that spends most of its 
time near the deterministic stable fixed points ±1. The right panel shows the stationary density for this 
process which has two distinct peaks. 

The interaction of the additive noise with the dynamics away from the deterministic limit cycle is 
the key to understanding the different switching behaviors. The two planar limit cycle systems we 
analyze have areas of the vector field where the rotation occurs in the opposite direction of the limit 
cycle. However, the two systems differ in the placement of this counterrotating behavior. In one system, 
the counterrotation occurs on only one side of the limit cycle, whereas in the other system, the limit 
cycle is surrounded on both sides by a counterrotating vector field. We will show that these different 
vector fields lead to very different bistable switching behaviors when perturbed by white noise. In 
one case, the noise gets amplified by the vector field, and switching between metastable states occurs 
frequently. However, in the other case, the noise gets suppressed, and switching between metastable 
states occurs very rarely. Through the use of stochastic phase reduction methods [21,24], we show that 
the directionality of switching is controlled by the drift of the phase reduced system while the rate of 
switching, similar to the Kramers’ problem, is determined by the magnitude of the diffusivity relative 
to the drift. Importantly, this work illustrates that a wide range of bistable switching dynamics can be 
modelled by limit cycle systems with additive noise, and provides further evidence that noise can affect 
seemingly similar nonlinear systems in strikingly different ways. 

The paper is organized as follows. First we present the equations for the two planar limit cycle 
systems we analyze. Next, we give the details of the transformation of the system to asymptotic phase 
which then allows us to reduce the two dimensional systems to a single stochastic differential equation 
for the evolution of the phase variable. We then review how the onset of bistable switching behavior 
can be predicted from the phase model by examining the drift term in the conservation (Fickian) form 
of the stationary Fokker-Planck equation. We then illustrate that even though the two systems display 
qualitatively similar stationary densities, the dynamics of the bistable switching in each system is quite 
different. By deriving asymptotic approximations to the stationary flux, we show that the directionality 
of the switching is determined by the sign of the flux. Lastly, by exploring the terms in the reduced 
phase system, we show that the rate of switching is determined by the relative magnitude of the drift 
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and diffusivity of the stochastic process. 


2 Models 

To explore this bistable switching behavior in limit cycle sytems, we consider the following analytically 
tractable deterministic oscillator 


x = -wy + 7*(1 - p 2 ) + cjyQ(p) = F(x,y) (2.1) 

y = wx + 72/(1 - p 2 ) - c~fxQ(p) = G(x,y), 

where p = x 2 + y 2 and the function Q(p) is such that Q(l) = 0 and determines the rotation away 
from the limit cycle. The dynamics of the system are easier to understand if (2.1) is rewritten in polar 
coordinates 


0 = w - "fcQ(p) (2.2) 

P = -1P{P 2 - !)■ 

In particular, it can now easily be seen how Q{p) affects the rotation away from the limit cycle at p = 1. 
We assume that the limit cycle is strongly attracting so that 7 > u is a large parameter. If c = 0, 
then (2.2) is independent of Q and yields the radial isochron clock model when 7 = 1 . In this case, the 
rotation 6 is constant away from the limit cycle and independent of p. If c / 0, then the rotation changes 
direction for values of p = p* such that Q(p *) = w/(c"f). We consider the following two cases: 

Q 1 (p)=p 2 - 1 (2.3) 

Q 2 (p) = w(1 -p) 2 . 

Qi yields the Stuart-Landau model [19]. However, for reasons made clear below, we will refer to the 
system with Q\ as the counterrotating (CR) whereas we will refer to Q 2 as the antirotating (AR) model. 
It follows that the limit cycle rotates counter clockwise with period 2tt/uj on p = 1 . 

Figure 2 plots the vector field of (2.1) for the two different choices of the function Q(p) (2.3). Note 
that, regardless of the definition of the function Q(p), the limit cycle is given by x(t) = cos(o jt) and 
y(t) = sin(wt) (thick black line). For the CR system, there is a unique value of p* = y/l + ^/(cy) where 
the rotation changes sign (dashed line); if c > 0 then p* > 1 and if c < 0 then p* < 1 (not shown). On 
the other hand, for the AR system, we have that p* = 1 ± yjljjcy). Hence, there are two values of p* 
if c > 0 (one inside the limit cycle and one outside, dashed lines) and none if c < 0. For simplicity, we 
only consider the c > 0 case for both systems. Note that for both the CR and AR systems, increases in 
c and 7 cause p* to move closer to the limit cycle. 

Next, we add independent white noise to each component 

x = F(x, y) + <JW x £ x (t) (2.4) 

V = G(x,y)+awy€ v (t), 

where a is the noise strength, w x and w y are constants, and £*(£), i = x,y, is Gaussian white noise with 
zero mean and unit variance (i.e., (£i(t)) = 0 and (£i(t)£i(s)} = S(t — s)). If w x = w y , then the stochastic 
process is rotationally symmetric. In this case, the CR system acts to amplify noise fluctuations whereas 
in the AR system, depending on the value of c, the noise can cause the system to rotate in the opposite 
direction of the deterministic limit cycle or even cease oscillating altogether [15]. When the statistical 
rotational symmetry is broken, w x 7 ^ w y , it was also shown that the CR and AR systems can display 
bistable switching between two metastable states ( [15] and see below). Since we are interested in further 
examining this bistable behavior, for the rest of the paper we assume, without loss of generality, that 
w x = 0 and w y = 1. This simplifies (2.4) to the following system of stochastic differential equations 
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Figure 2: Vector fields of the two planar systems. Thin lines represent streamlines of the deter¬ 
ministic vector fields. Thick line is the limit cycle, (a) The counterrotating (CR) system with c = 4 and 
7 = 1. The vector field changes rotation direction at one value of p outside of the limit cycle (dashed 
line), (b) The antirotating (AR) system with c = 15 and 7 = 5. The vector field changes rotation 
direction at two values of p one outside and one inside of the limit cycle (dashed lines). 


(SDEs) 


x = F(x,y) (2.5) 

V = G(x, y) + cr£(t), 

To simulate the SDEs, we employ the Euler-Mayurama method with a sufficiently small time step. 

3 Bistable Switching in Limit Cycle Systems 

In previous work, it was shown that with c ~ 1 and moderate levels of noise ex, the two planar limit cycle 
systems exhibit bistable-like switching. Figure 3 plots example trajectories for both systems. For the CR 
system, the x coordinate appears to display strong bistable switching between ±1 at a fairly frequent 
rate. However, for the AR system, it is the y coordinate that appears most like a bistable switch, but 
does so at a rate much slower than that of the CR system. This shows that the vector field away from 
the limit cycle interacts with the noise in interesting ways in order to produce this behavior. Thus, we 
seek to understand (i) how this bistable behavior emerges in each of these systems, and (ii) why bistable 
switching is different between the two systems. To do so, we first utilize a phase reduction technique for 
stochastic limit cycle systems [21,24] in order to reduce the system to a one-dimensional Ito SDE which 
describes the position of the system on its underlying determinsitic limit cycle. 

3.1 Reduction to a Phase Model 

In the absence of noise both planar sytems (2.1) exhibit stable limit cycle oscillations. Since we also 
assume that the rate of attraction back to the limit cycle is large, i.e., 7 w, this places us in the 
ideal situation use so-called phase reduction techniques [6, 13,21,24]. In this case, the dynamics of a 
single phase variable ip that denotes the location of a trajectory on the underlying limit cycle will be an 
accurate (to order I/ 7 ) description of the full dynamics of the system. 
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Figure 3: Bistable switching in the CR and AR systems. Dashed lines show the deterministic 
limit cycle a = 0 (which is the same in both cases). For each system 7 = 50 and ui = 1 (a) CR system 
with c = 4 and a = 0.5. Bistable switching is most prominent in the x-component. (b) AR system with 
c = 20 and a = 0.45. Bistable switching is most prominent in the y-component. 


Consider the noiseless system (2.1) which has a stable 27r/w-periodic solution. We define a phase 
<p £ [ 0 , 2tt/uj) coordinate on the limit cycle solution such that it increases by 2it/u> for every cycle of the 
system. That is, p is defined by 


dp = dpdx dpdy = 

dt dx dt dy dt 

We use phase instead of polar coordinates because, in the absence of external input, the phase variable 
ip increases monotonically with a constant rate. Thus, any change to the rate of increase of the phase 
variable will solely be due to the external input (in this case noise) and not due to the intrinsic dynamics 
of the system. This is not the case for the angle variable 0 in polar coordinates (2.2). 

In general, it is not possible to find the exact mapping from Cartesian coordinates to <p. However, 
for the planar systems we use here, the mapping ( x , y) 1 —>■ ip can be found analytically. To see this, recall 
the polar coordinate form of the systems (2.2). Owing to the simplicity of this system, we can easily 
transform the angle 6 to tp using 


1 r 


<P = 


6 + 7 cQ{p) 


Integrating both sides of the above equation with respect to time yields 




1 

w 


6 + 7 c 


1 

w 


9 + 7 c 


Q(p)dt 


Q(p ) 

-IPiP 2 - 



Note that we have ignored the integration constant which is an arbitrary phase shift related to the initial 
starting point. Next, we define a new amplitude variable r = p — 1 which is the distance of the system 
from the deterministic limit cycle. Using the fact that 9 = tan ~ 1 (y/x) and p = \J x 2 + y 2 , we arrive at 
our desired mapping 


ip = — [tan l (y/x) + off(r)] , r = \Jx 2 + y 2 — 1 , 


(3.2) 
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where H(r) = f _ ( - T . + ^^ 7 ^ L jp_ 1 - ) dr. It is important to note that the function H(-) depends on the 
function Q(-). Thus, the mapping to phase will be different for the CR and AR systems. 

Lines of constant <p are called isochrons [23] - all deterministic trajectories starting on an isochron 
converge as t —» oo. Hence, isochrons encode information about how the phase of a trajectory responds 
to a perturbation away from the limit cycle. Indeed, the gradient of this mapping i = x, y is typically 
referred to as the phase response curve of the i-th component of the limit cycle [6, 13, 18] as it quantifies 
the change in phase due to a d-function perturbation at a particular phase on the limit cycle. 

With our mapping in hand (3.2), we transform the stochastic Cartesian system (2.5) to (tp,r) coor¬ 
dinates 


ip = 1 + yn(i/;,r) + ah(<p,r)£(t) (3.3) 

(^2 

r = 7/(r) + —n r {y,r) + ag(ip,r)(,(t), 

where, 


h &> r ) = ^ = uirW) Sin(a + ^ ’ 9 ^ r) -% = Sin(a) 

dr dr 

f(r ) = T^_F(x(ip, r), y(ip, r)) + —G(x(ip, r),y(ip, r)) = -r(r + l)(r + 2) 

A = c(r + l)H'(r), %[} = tan _ 1 (A), a = uiip + cH (r), 

and we have used the inverse mappings x((p, r) = (r + 1) cos(a) and y(ip, r) = (r + 1) sin(a). While the 
noise is additive in cartesian coordinates (2.5), it is multiplicative in (ip,r) coordinates. As such, we use 
the Ito interpretation. The terms that arise from the stochastic change of variables (i.e., the terms that 
come from the Ito calculus [ 8 ]) are n(ip,r) = h(<p,r)^ + and n r ((p,r) = h(ip,r)^ + 

We stress that (3.3) is an exact transformation, and no approximations have been used up to this point. 
Thus, our goal now remains to project out the variable r in order to arrive at a single differential equation 
for the phase. 

Assuming that the deterministic limit cycle is highly attracting (i.e., 7 cu), the system should cling 
tighly to the underlying deterministic limit cycle, and on average r will be approximately 0. Thus, we 
can set r = 0 in (3.3) to obtain (to leading order in 7 -1 ) the one dimensional approximation 

2 

~ 1 + 7 »w, (3.4) 

where no(ip) = n(ip, 0 ) and ho(tp) = h(ip, 0 ) is the phase response curve of the y-component of the 
deterministic limit cycle. Note that (3.4) can be rigorously derived by using adiabatic (or quasi-steady 
state) reduction techniques [ 8 , 21 ], 

The stationary density u ss (ip ) of the reduced system (3.4) will provide insight into the dynamics of 
the full model. The stationary density obeys the following Fokkcr-Planck (FP) equation 


0 — dip ^dip\^Dip{ip^u ss ^ L^((p)'u ss ], (3.5) 

with periodic boundary conditions. Note that in general analytical solutions to the above FP equation 
are not possible and we instead find approximate solutions using a forward or backward finite-difference 
scheme (depending on parameters). = 1 + ^-no{ip) is the Ito drift, or the deterministic portion of 

the SDE (3.4), and Dip(ip) = is the diffusivity. The FP equation (3.5) can also be written in 

conservation or Fickian form 


0 — dp [Dip^tp'ydpUss fjp((p)'u ss ], 

where V v (ip) = I v {ip) — D[’(<p) is the Fickian drift. 


(3.6) 
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The Fickian drift is closely related to the stationary density of the the system. To see this, integrate 
both sides of (3.6) to obtain 




where J is the probability flux. If the flux is zero, then u ss 


has the solution u ss (<p) = Ne where 


$(A) 



-v(<p) 

D {<f>) 


dip, 


and N is a normalization factor. 

On the other hand, the Ito drift represents the average increment of the stochastic process during an 
infinitesimal time dt. For a periodic system such as ours, the Ito drift is also related to the stationary 
probability flux. To see this, integrate both sides of (3.5) to obtain 


J — dip (Dip{jp)llss') I(p(}P')V J SS, 

where J is again the probability flux. Lastly, integrating both sides of the above equation from 0 to 27r/w 
(and noting that D v (<p) and u ss are 2tt/uj periodic) yields 

l 1 r 2 -k/ui 

J = - 2/kJuj J 0 I v(v) u ss(ip) d ip- (3-7) 

In the case that the stationary density is uniform u ss = (27r/w)^ 1 , the probability flux is proportional to 
the average Ito drift over one period of the oscillations. However, in general the stationary density will 
not be uniform, and its shape will have a large influence on the flux. 

The flux can also be thought of as the average frequency of the oscillations. In the noiseless system, 
Ip(p) = 1, and thus the flux is (27 t/oj) _ 1 , which is the frequency of the deterministic limit cycle. In fact, 
for stochastically perturbed periodic systems that have an Ito drift of 1, the stationary flux is always 
(27r/o;) _1 regardless of the noise level. 


3.2 Diagnosing Bistability in the CR and AR Systems 

Before showing how the onset of bistable switching behavior can be predicted using the phase reduced 
system, we first give the details of the drift and diffusivity functions which will be vital to our later 
analysis. For the CR system, the Ito drift, Fickian drift, and diffusivity are given by 


lfp R {p) = 1 — [ccos(2w<p) + sin(2o;</j)] (3.8) 

2 

V£ R (ip) = 1 + [—c 2 sin(2w</?) + ccos(2w(/?)] 

2 

Dp R (<p) = [c 2 sin 2 (ojp) — csin(2wv3) + cos 2 (uj<p)]. 

The corresponding quantities for the AR system are 


t ar 


M = 1 ~ ^r 


- sin 2 (w^) d-sin(2w</?) 

2 uj 


V2 R (P) = 1 - —csin 2 (unp) 


Dp R {p) = — cos 2 (uip). 


(3.9) 
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Figure 4 plots the above quantities along with the stationary density of the phase reduced system (nu¬ 
merical solution of equation (3.6)) and the stationary density computed from Monte-Carlo simulations 
of the full system (3.3). Note the excellent agreement between the phase reduction and the full system 
(compare open circle to the solid lines in the lower panels). The CR (AR) system is shown in Fig. 4 
(a) and (c) ((b) and (d)) with two different values for c and all other parameters held constant. Recall 
from section 2 that increases in the parameter c causes the oppositely-rotating region(s) of the vector 
field to move closer to the limit cycle (i.e., increasing c brings p* closer to 1 ). It can clearly be seen 
that, in both systems, increasing the value of c also causes the Fickian drift to become negative which 
signals the onset of the bistable switching behavior. Furthermore, the two points where V v changes sign 
from positive to negative (see upper panels of Fig. 4 (c) and (d)) as p is increased are local maxima of 
the stationary density, where the oscillator spends most of its time, struggling to move past a dynamical 
barrier imposed by the noise. Thus, the zeros of the Fickian drift that have negative slope are metastable 
points that the system switches between. It is important to note that the noise-induced contributions to 
the Fickian drift (the a 2 terms) are what causes it to become negative. Indeed, if a is close to zero, then 
V v ss 1 for both systems. 


CR System 

(a) 


AR System 

(b) 




Figure 4: Bistability occurs when the Fickian drift becomes negative. Upper panels plot the 
Fickian drift (V^, solid lines), Ito drift (J^, grey dash-dotted lines), and diffusivity (D v , solid lines) for 
the phase reduced system (3.4). Lower panels plot the stationay density of the phase reduced system 
(solid lines) obtained from solving (3.5) along with Monte-Carlo simulations of the full system (3.3) 
(open circles) with 7 = 50. In all cases, u> = 1. The (a) CR system with c = 2 and <7 = 0.5 and the 
(b) AR system with c = 10 and a — 0.55 both have relatively flat stationary densities. Strong peaks in 
the stationary densities occur when the Fickian drift crosses zero with a negative slope in (c) CR system 
with c = 4 and er = 0.5 and (d) CR system with c = 20 and a = 0.5 


























We can already see interesting differences between the two systems emerge in their respective phase 
reductions. In the AR system, the diffusivity is independent of c which is not the case for the CR system. 
Thus, in the CR system, changing the dynamics of the vector field around the limit cycle (by varying 
c) affects both the drift and diffusivity of the phase reduced system, whereas in the AR system, the 
fluctuations of the phase variable are unaffected by the dynamics off the limit cycle (compare the dashed 
traces). It can also be seen that in the AR system, the diffusivity is small compared to the magnitude of 
the drift (compare the dashed lines to the solid lines in the upper panels of Fig. 4 (b) and (d)). Lastly, 
in the AR system, the the Ito and Fickian drifts are approximately the same whereas in the SL system, 
they are quite different (compare the dash dotted lines to the solid lines in the upper panels of Fig. 4). 

Thus, the noise interacts with the nonlinear dynamics of the CR and AR systems in different ways 
to produce stochastic phase variables with qualitatively similar stationary behavior. However, recalling 
the trajectories shown in Figure 3, and noting the differences in the shapes of the stationary densities in 
the lower panels of Figure 4 (c) and (d), we can begin to see that there may be stronger differences in 
the transient behavior of the CR and AR systems. Indeed, if we plot stochastic phase trajectories from 
simulations of the full system (3.3) as in Figure 5, we can see strong differences emerging between the 
two systems. The first thing to note is that the switching behavior observed in Cartesian coordinates 
translates to sudden jumps in the trajectory of p that we will refer to as “phase jumps”. The CR system 
displays fairly frequent phase jumps whereas in the AR system, these transitions are rare (compare the 
x-axis scale in Fig. 5 (a) and (b)). This is similar to what we observed in simulations of the Cartesian 
system in Figure 3. Secondly, in the CR system, the phase jumps are in the positive phase direction 
(the same direction as the deterministic limit cycle) but in the AR system the jumps are in the opposite 
direction of the deterministic limit cycle. 



Figure 5: CR and AR systems display different switching dynamics. Sample phase trajectories 
from simulations of the full system (3.3) with 7 = 50 and w = 1. (a) CR system with c = 4 and a = 0.5 
displays frequent switching while (b) AR system with c = 20 and a = 0.55 displays very rare switiching 
(note difference in x-axis scale). 

We now seek to understand why the bistable behavior differs between the two systems. We will 
first concentrate on the different switching directions. In particular, we will show that the stationary 
probabilty flux is related to the directionality of the rotation of the oscillators, and we will derive 
analytical approximations to the flux in order to gain insight into the mechanism underlying the different 
switching directions. 


4 Stationary Probability Flux Determines Switching Direction 

Here, we derive analytical approximations to the flux for the AR and CR systems. For both systems, 
the bistable behavior increases as c and a are varied such that the Fickian drifts becomes more negative. 
Because of this ambiguity, we introduce different scalings for the two systems in order to arrive at a 
single parameter which, when increased, causes bistable switching to occur. With this scaling, we are 
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also able to take a distiguished limit (c —> oo and a —»■ 0 ) of the phase reduced systems in order to obtain 
analytically tractable approximations. We note that if c becomes too large, then the phase reduction itself 
may cease being an accurate description of the full system. However, as long as we assume 7 > c > 1, 
then the results in the next section hold, and we verify this fact with the use of numerical simulations. 


4.1 CR System 

We first introduce the parameters £ = \ and Kcr = f ■ Keeping the parameter Kcr constant as we 
simultaneously let a —>■ 0 and e —> 0 allows us to expand the system in an asymptotic series in e. To see 
this, substitute c = \ and a = eKcr into equation (3.8) to obtain 


= 1 ~ cos(2ojip) — e 2 ^ CR sin(2u;y>) (4.1) 

v 2 2 uj 

EE ^(^+^1^)+^) 

V£ R {<P) = 1 - 4^ sin(2u^) + £^^ cos(2w<^) 

= V£ Rfi {ip)+eVg R +(ip) 

= ^ CR sin 2 (w</j) — £ K CR sin(2w</j) + £ 2 ^~ CR cos 2 (wy>) 
v 2 w y J 2 w ’ 2w v ’ 

= DC R ’°(<p)+eDC R ’\v) + 0(e 2 ). 


Next, we expand the stationary density u ss = u® s + ul s + 0(e 2 ) to obtain an asymptotic approximation 
to the flux (3.7) 


J = 


2tt/u 


(-2tt / oj 


1 + £ 


tCR, 1 


(vhss^dtp 


0(s 2 ), 


(4.2) 


where we used the fact that I^ R '°(ip) = 1, u° ss (<p)dtp = 1, and u] s (<p)dip = 0. Thus, in order 

to close our approximation, we must find the leading order approximation to the stationary density u® s . 

Starting from the Fickian form of the FP equation (3.6), we integrate once and expand terms in 
powers of e to obtain the leading order equation 


-J° = D^ R ’°( V )d v u° at - V? R ’ 0 {<p)u° aa , 

with u° s (0) = u® s (2tt/u]) and J° = ( 2tt/oj)~ 1 is the leading order approximation to the flux. Owing to the 
simple forms of V^ R ’° and D^ R, °, we integrate the above equation and apply the boundary conditions 
to obtain 


zt° s = ^ e [e Al Ei(— Acot(u<p) + Ai) — e Al Ei(—Hcot(w(/?) — Ai)\ , (4.3) 

27r/o; J^cr 

where A = 2u>/Kq R , 

r _yCR,o f£\ 

= / Jjio,-' d( P = M sin V^)) + Acot{wy), 

J Dp ’ (ip) 

and Ei(x) is the exponential integral. Figure 6 plots approximation u° ss (4.3) against Monte-Carlo 
simulations of the full model (3.3) for three different values of Kcr■ Notice that bistability increases 
as Kcr is increased and that our approximation is in excellent agreement with the simulations for all 
values of Kcr■ Thus, bistability occurs when the product of ca is large enough, illustrating that the noise 
interacts with the vector field in a nonlinear way in order to produce the bistable switching behavior. 

Now that we have an approximation to the stationary density, we can compute our approximation to 
the stationary flux (4.2). Unfortunately, the integral in equation (4.2) can not be computed analytically, 
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Figure 6: Approximations of the CR stationary density. The asymptotic approximation (4.3) 
(solid line) is plotted against Monte-Carlo simulations of the full system (3.3) (open circles) for (a) 
Kqr = 1, (b) Kcr = 2 , and (c) Kcr = 4. In all figures, 7 = 50, c = 20, and oj = 1 . 


and we therefore approximate it using quadrature. Figure 7 plots the approximation to the flux (solid 
line) along with the results from Monte-Carlo simulations of the full system (3.3). The dashed line is 
( 2tt/uj )~ 1 which is the leading order approximation to the flux. The lower panels plot example phase 
trajectories for three different values of Kcr- To compute the flux numerically, we took advantage of 
the fact that the system is ergodic and simulated (3.3) for a long time T. Using the time series ip(t) 
obtained from the simulation, we then computed the time average of the Ito drift 


1 1 


2'k/loT J 0 


It can clearly be seen that the simulations are in excellent agreement with our approximation. As Kcr 
is increased, the flux decreases and eventually becomes negative. The value of Kcr where the flux 
becomes negative is also the value of Kcr where the Ito drift starts to become negative. Ignoring the 
e 2 terms in I^ R (4.1), we find that the value of Kcr where the Ito drift first becomes negative is given 
by 1 Kcr = \J/e which is approximately 6.32 for the parameters in Figure 7. Thus, we can see that 
around this point, the flux is close to zero and eventually becomes negative as Kcr is increased further. 

When the flux is positive, the phase jumps in the direction of the deterministic limit cycle (see bottom 
left panel in Figure 7). This explains the behavior we saw in Figure 5 (a) as Kcr = 2 for the parameters 
used. However, when the flux is close to zero, the phase still jumps fairly regularly, however the jumps 
appear to be equally likely in both the postive and negative directions (bottom middle panel in Figure 
7). Lastly, when the flux becomes negative, the process jumps in the direction opposite of the motion of 
the deterministic limit cycle (bottom right panel in Figure 7). Thus, it is clear that for the CR system, 
the flux determines the average direction of the phase jumps, but does not influence the rate at which 
the process jumps. 


4.2 AR System 

As in the previous section, we introduce the parameters e = 4 and Kar = a e ■ Note the difference 
between Kcr and Kar- Substituting c = t and a 2 = sKar into equation (3.9), we obtain 

1 Set min^, I^ R ’°(ip) + Ip R,1 (<p) = 1 — e K ^J i = 0 and solve for Kcr- 
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Figure 7: Approximation of the CR flux. The asymptotic approximation (4.2) (solid line) is plotted 
as a function of Kcr along with simulations of the full system (3.3) open circles. The dashed line is 
(27r/w) -1 which is the frequency of the deterministic limit cycle and the zeroeth order approximation to 
the flux. Lower figures are sample phase trajectories of the full system (3.3). From left to right Kcr = 4, 
Kcr = 6 . 6 , and Kcr = 7.5. In all figures, 7 = 50, c = 20, and oj = 1. 


I$*{<P) = l-^sin 2 (^)- £ ^sin( 2 ^) (4.4) 

V* R (<p) = 1 - ^^sin 2 (w<p) 

D^ R {y) = £ ~^J~ cos 2 (w</?). 

With this scaling, it is very easy to see that l£ R « V£ R and that D^ R ss 0. Using equation (3.6), 
and noting that the leading order diffusivity is 0 , we obtain a very simple equation for the leading order 
stationary density u® s 


u° ss (f) 


Jo 

vJWY 


(4.5) 


From the above equation, we can see that the leading order flux is the same as the normalization constant 
for the density. Enforcing that u® s ((p)d(p = 1 we find 


_ a /4 — Kar 
A n ju 1 


(4.6) 
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However, notice that this approximation is only valid for Kar < 4. This is an interesting point that we 
will return to in the next section. Figure 8 plots our approximation u® s against Monte-Carlo simulations 
of the full model (3.3) for three different values of Kar- Again, notice that bistability increases as Kar 
is increased. However, notice that our approximation begins to break down near Kar = 4 (Fig. 8 (c)), 
but this can be alleviated by taking 7 to be an order of magnitude larger (compare circles to crosses). 
This change in the behavior of the system can be seen even more clearly when we plot the flux as a 
function of Kar as in Figure 8 . In contrast to the CR system, the flux in the AR system quickly goes 
to zero as Kar is increased. We can see that the analytical approximation works very well for values of 
Kar away from 4 (compare circles to solid line). Taken together with Figure 8 (c), this shows that, near 
Kar = 4, the 0{ I/ 7 ) terms that we are ignoring in the phase reduction (3.4) become significant. After 
Kar = 4, our approximation (4.5) is no longer valid. However, numerical simulations show that the flux 
is very close to zero. 


(a) 


(b) 


( c ) 





Figure 8 : Approximations of the AR stationary density. The asymptotic approximation (4.5) 
(solid line) is plotted against Monte-Carlo simulations of the full system (3.3) (open circles) for (a) 
Kqr = 2, (b) Kcr = 3.5, and (c) Kqr = 3.9. In all figures, 7 = 50, c = 20, and w = 1. The crosses in 
(c) are simulations of the full system with 7 = 50. 

As with the CR system, we can see that the flux is a good indicator of the average direction of the 
phase jumps (see the bottom panels in Fig. 9). However, it is difficult to interpret the results in the 
Kar > 4 regime as the flux is very close to zero. In contrast to the CR system, the rate of switching 
in the AR system also decreases with Kar- Indeed, at Kar = 3.5, the system displays frequent phase 
jumps in the positive direction, but at Kar = 6 the system displays incredibly rare jumps in the negative 
direction (compare x-axis scales in lower left and lower right panels) which is similar to the behavior 
we saw in Figure 5 (b) as Kar = 6 for the parameters used. However, this decrease in the jump rate 
cannot be due to the flux, as we saw that with near zero flux in the CR system, the rate of switching 
was still fairly frequent. In the next section, we seek to explain why there are differences in the switching 
frequency. 


5 Switching Frequency is Determined by Diffusivity 

We have shown that the onset of bistable switching in the CR and AR systems occurs when the Fickian 
drift becomes negative, and that the direction of the switching (either in the same or opposite direction 
of the deterministic limit cycle) is indicated by the sign of the stationary probability flux. Here we will 
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Figure 9: Approximation of the AR flux. The asymptotic approximation (4.6) (solid line) is plotted 
as a function of Kar along with simulations of the full system (3.3) open circles. The dashed line 
is (27t/o;) _1 which is the frequency of the deterministic limit cycle. Lower figures are sample phase 
trajectories of the full system (3.3). From left to right Kar = 3.5, Kar = 3.9, Kar = 4.5, and 
Kar = 6 . In all figures, 7 = 50, c = 20 , and oj = 1 . 


show that the frequency of the switching events can be explained by the relative strengths of the Ito 
drifts and diffusivities in the reduced phase system (3.4). 

5.1 Rare Switching in AR System is Caused by Weak Diffusivity 

To understand both the decreasing rate of switching and the change in switching direction of the AR 
system, we examine the contributions of the Ito drift and diffusivity to the reduced phase model using 
the scalings introduced in the previous section. First, recall the diffusivity for the AR model given in 
(4.4). Notice two things about the AR diffusivity: (1) that it is order e; and (2) that it has two zeros 
ip* = A, regardless of the values of Kar and e. We will refer to these zeros of the diffusivity as 
“ratchet points”. The only way that the reduced phase system can cross these ratchet point is if the Ito 
drift is nonzero at these points. Figure 10 (a) plots the Ito drift (solid) and diffusivity (dashed) as a 
function of ip in the left panel and a sample stochastic phase trajectory of the full system (3.3) in the 
right panel when Kar = 3.5. It is clear that the Ito drift is always significantly larger than the diffusivity, 
and thus is the main determinant of the motion of the stochastic process. This can be seen in the right 
panel of Fig. 10 (a) where the phase variable increases steadily over time. The effect of the changing 
Ito drift can be seen in the slope of the trajectory as the phase variable increases more slowly when it is 
near the two ratchet points. In this case, the Ito drift is positive at the ratchet points (as indicated by 
the horizontal arrows) which indicates that the stochastic trajectory will pass over the ratchet points in 
the direction of the deterministic limit cycle. 

As Kar is increased to 3.961 as in Fig. 10 (b), the Ito drift near the ratchet points becomes the same 
order at the diffusivity (inset, left panel). Thus, the process begins to have a more difficult time passing 
over the ratchet points, and begins to display the phase jumps seen previously (right panel). The jumps 
occur in the positive direction as the Ito drift at the ratchet point is still positive (horizontal arrows, left 
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panel). When Kar = 4 (Fig. 10 (c)), an interesting phenomenon occurs: the Ito drift is exactly zero at 
the ratchet points (to see this, plug ip* = into the equation for in (4.4)). Thus, the phase 

reduced system (3.4) predicts that once a trajectory arrives at a ratchet point, it will be stuck there 
indefinitely. However, it is clear from the 7 = 50 sample trajectory in the right panel of Fig. 10 (c) that 
the 0 ( 1 / 7 ) terms in the phase reduction become significant as the trajectory of the full system (3.3) is 
able to pass over the ratchet points. This discrepency between the full system and the phase reduced 
system is caused by the fact that the phase reduction ignores amplitude fluctuations. That is, for these 
values of Kar and 7 , fluctuations in the amplitude variable allow the full system (3.3) to move past 
the ratchet points. However, as 7 is increased, the trajectory has a more difficult time overcoming the 
ratchet points (7 = 500 and 7 = 5000 traces). Thus, as 7 —> 00 the stationary density for ip converges to 
a delta function at one of the two ratchet points depending on the initial conditions. 

Lastly, when Kar > 4 (Fig. 10 (d)), the Ito drift becomes negative at the ratchet points, and the 
phase reduction predicts that the process can only move past the ratchet points in the opposite direction 
of the deterministic limit cycle. However, the Ito drift is only negative in a neighborhood of the ratchet 
points, and is strongly positive everywhere else. Furthermore, the diffusivity is still much smaller than 
the Ito drift. Thus, away from the ratchet points, the Ito drift dictates that the process move in the 
same direction of the limit cycle. Therefore, the phase reduction predicts that the only way that the 
process can jump between ratchet points is via a rare event where diffusion overcomes the positive Ito 
drift and the process jumps in the opposite direction of the deterministic limit cycle (right panel, 7 = 50 
trace). Thus, the rare switching in the AR system is caused by the fact that the Ito drift is negative near 
the ratchet points, but positive and larger than the diffusivity away from the ratchet points. In contrast 
to the phase reduced system, if 7 is made smaller, fluctuations in the amplitude variable allow the full 
system (3.3) to pass over the ratchet points in the same direction of the deterministic limit cycle (right 
panel, 7 = 10 trace). Note also that the small diffusivity is the reason why our asymptotic approximation 

(4.5) breaks down when Kar > 4. That is, our approximation assumes that the diffusivity is 0 for all <p. 
However, close to the points where the Ito drift first becomes negative, the diffusivity is non-negligible 
(and is in fact the reason why a switch can actually occur). 

We can quantify the accuracy of the phase reduced system (3.4) by computing the L 2 norm of the 
difference between the stationary density obtained from solving the FP equation for the reduced system 

(3.5) u ss and the stationary density obtained from Monte-Carlo simulation u^ c . More precisely, we take 
the error to be 


Error = 



(5.1) 


Figure 11 plots this error as a function of Kar and two different values for 7 . Note that to solve (3.5), 
we found it best to use a forward finite difference scheme when the flux was positive, and a backwards 
finite difference scheme when the flux was negative. It is clear that the phase reduced model is highly 
accurate for Kar < 4, but the error quickly increases as Kar —► 4+. Interestingly, the 0 ( 1 / 7 ) terms 
appear to be more important for accuracy when Kar > 4. Regardless, increasing 7 causes the error to 
decrease for all values Kar 7 ^ 4 (compare solid and dashed lines). 


5.2 Faster Switching in the CR System is Caused by Diffusivity Being Larger 
than the Ito Drift 

For the CR system, the ratchet points (zeros of the diffusivity) occur at approximately p* = 0 , 7 r as 
Dp R (tp*) = 0(e 2 ) (see (4.1)) independent of Kqr • The left panel Figure 12 (a) plots the Ito drift and 
diffusivity for the CR system when Kqr = 4. In this case, the diffusivity is much larger than the Ito drift 
except near the two ratchet points. Thus, switching between the metastable states (zeros of the Fickian 
drift with negative slope) will occur at a higher rate than in the AR system. At the ratchet points, the 
Ito drift is positive, which indicates that the reduced phase system (3.4) predicts that the process moves 
past the ratchet points in the direction of the deterministic limit cycle. A sample trajectory from the 
full system (3.3) in the right panel shows that this is indeed the case as ip increases with time. 
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Figure 10: AR system is drift dominated. Left panels show the Ito drift (solid lines) and diffusivity 
(dashed lines) of the phase reduced system (3.4). Solid vertical lines indicate the ratchet points (where 
the diffusivity is zero) and the horizontal arrows indicate the sign of the Ito drift at the ratchet points. 
Right panels show sample phase trajectories of the full system (3.3). For all plots, c = 20, u> = 1, and 
unless indicated otherwise 7 = 50. (a) Kar = 3.5, (b) Kar = 3.961, (c) Kar = 4, and (d) Kar = 6 . 

Figure 12 (b) shows the behavior of the CR system when Kqr = \f2uj/e ss 6.32 where = 

0(e 2 ) (see (4.1)). Thus, the Ito drift is approximately zero at the ratchet points. As in the AR system, the 
phase reduction predicts that once a trajectory arrives at a ratchet point, it will be stuck there indefinitely. 
However, the sample trajectories in the right panel again show that accuracy of this predicition is highly 
dependent on the value of 7 . Lastly, when Kar > y/^uj/e (Fig. 10 (c)), the Ito drift becomes negative 
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Figure 11: Error in the phase model approximation of the stationary density for the AR 
system as a function of Kar ■ L 2 norm of the difference between the stationary density obtained from 
solving (3.5) and the stationary density obtained from Monte-Carlo simulations of the full system (3.3) 
as a function of Kar when c = 20, u> = 1, 7 = 50 (solid line), and 7 = 500 (dashed line). 


at the ratchet points, and the phase reduction predicts that the process can only move past the ratchet 
points in the opposite direction of the deterministic limit cycle. In contrast to the AR system, the CR 
system has a diffusivity that is larger (in magnitude) than the Ito drift away from the ratchet points. 
Thus, the switching in the opposite direction on the limit cycle will occur at a higher rate than in the 
AR system (compare right panel of Fig 12 (c) to right panel in 10 (d)). Just as in the AR system, if 
7 is made smaller, the full system can also overcome the ratchet points in the same direction of the 
deterministic limit cycle by moving off the limit cycle (right panel, 7 = 10 trace). 

Just as we did with the AR system, we quantify the accuracy of the phase reduced system (3.4) by 
exploring the error (5.1) between the stationary density obtained from solving the FP equation for the 
reduced system (3.5) u ss and the stationary density obtained from Monte-Carlo simulation u^l c . Figure 
13 plots this error as a function of Kcr and two different values for 7 . Interestingly, the error overall 
appears to be larger than what we found in the AR system (compare y-axis scale in Figures 13 and 11). 
Nonetheless, the qualitative trend is the same in both systems: as Kcr approaches the value where the 
Ito drift is zero at the ratchet points, the error in the phase reduction increases. Furthermore, increasing 
7 causes the error to decrease for all values Kcr 7 ^ a/ 2 w/e (compare solid and dashed lines). 

6 Discussion 

We have shown that a limit cycle system with additive noise can display a wide range of bistable switching 
dynamics. Using stochastic phase reduction methods [21,24], we are able to relate different aspects of 
this bistable switching behavior to terms appearing in the phase reduced model. That is, the emergence 
of bistable switching in a limit cycle system occurs when the drift term of the conservation (Fickian) form 
of the Fokkcr-Planck equation for the phase reduced system becomes negative. The zeros of the Fickian 
drift with negative slope represent metastable points that the stochastic process switches between and 
correspond to peaks in the stationary density for the phase. The rate and directionality of the switching 
are controlled by the drift term in the stochastic differential equation for the phase reduced system, i.e., 
the Ito drift. We have shown that the Ito drift determines the stationary probability flux, and thereby 
the directionality of the switching. That is, if the flux is positive (negative), the system will (on average) 
switch between the metastable points in the same (opposite) direction of the deterministic limit cycle. 
Lastly, as with the Kramer’s problem, the magnitude of the Ito drift relative to the diffusivity determines 
the rate at which the process switches between the metastable points. 
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Figure 12: CR system is diffusion dominated. Left panels show the Ito drift (solid lines) and 
diffusivity (dashed lines) of the phase reduced system (3.4). Solid vertical lines indicate the ratchet 
points (where the diffusivity is zero) and the horizontal arrows indicate the sign of the Ito drift at the 
ratchet points. Right panels show sample phase trajectories of the full system (3.3). For all plots, c = 20, 
u> = 1, and unless indicated otherwise 7 = 50. (a) Kqr = 4, (b) Kcr = 6.32, and (c) Kcr = 7.5. 

We illustrated these findings with the use of simple planar limit cycle systems that have areas in their 
vector fields where rotation occurs in the opposite direction of the limit cycle. In the CR system, this 
counterrotation occurs on only one side of the limit cycle. This lead to the diffusivity the phase reduced 
system being much larger in magnitude than the Ito drift. Thus, the CR system displays fairly frequent 
swithching between metastable points. In the AR system, the counterrotation occurs both inside and 
outside of the limit cycle. The phase reduction for the AR system had a diffusivity that was much smaller 
in magnitude than the Ito drift, and thus switching between metastable points is a rare event. Both 
systems can display switching in either the same or opposite direction of the deterministic limit cycle. 
Thus, it is the differences in the dynamics off the limit cycle between the two systems that is responsible 
for the different rates of switching. This highlights the fact that in stochastic limit cycle systems, if the 
interactions of the noise with the vector field off the limit cycle are ignored in the phase reduction (i.e., 
the 0(a 2 ) terms in (3.4)), then the phase model will fail to capture many interesting phenomena that 
are present in the full system ( 2 . 1 ). 

We have also identified an interesting singularity where the leading order phase reduced system fails 
to capture the dynamics of the full system regardless of the rate of attraction back to the limit cycle. 
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Figure 13: Error in the phase model approximation of the stationary density for the CR 
system as a function of Kcr■ L 2 norm of the difference between the stationary density obtained from 
solving (3.5) and the stationary density obtained from Monte-Carlo simulations of the full system (3.3) 
as a function of Kcr when c = 20, u> = 1, 7 = 50 (solid line), and 7 = 500 (dashed line) 


This occurs when the Ito drift of the phase reduced system has a zero at exactly the same point as the 
diffusivity. In this case, the phase reduction predicts that the system should remain stuck at this point 
for all time. However, simulations of the full system reveal that this is not the case. Thus, 0 ( 1 / 7 ) terms 
that are ignored in the leading order reduction become necessary to capture the behavior of the full 
system. 

This work identifies an alternative mechanism for the emergence of bistable switching dynamics 
observed in a number of physical systems ranging from the dynamics of genetic swiches in bacteria [9,14, 
22] to up and down states in the cortex thought to be involved in working memory [5,7]. Interestingly, 
the dynamics of the phase variable in the CR system (Figure 5 (a)) is reminiscent of the behavior of 
a brownian ratchet, which is a standard model for molecular motors [1,12,20]. Furthermore, as this 
bistable switching behavior has also been shown to occur in more complex limit cycle systems [15], 
this could lead to interesting consequences for coupled oscillator systems. Indeed, an interesting type 
of mixed synchronization has already been observed when oppositely rotating oscillators are coupled 
together (see, for example, [3,17]). Exploring the behavior of coupled bistable oscillators could lead to 
new synchronization or phase-locking phenomena. 
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